Vortex crystallisation in classical field theory 
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We show that the formation of a vortex lattice in a weakly interacting Bose condensed gas can be 
modeled with the nonlinear Schrodinger equation for both T = and finite temperatures without the 
need for an explicit damping term. Applying a weak rotating anisotropic harmonic potential we find 
numerically that the turbulent dynamics of the field produces an effective dissipation of the vortex 
motion and leads to the formation of a lattice. For T = this turbulent dynamics is triggered by 
an already known rotational dynamic instability of the condensate. For finite temperatures, noise is 
present at the start of the simulation and allows the formation of a vortex lattice at a lower rotation 
frequency, the Landau frequency. These two regimes have different vortex dynamics. We show that 
the multimode interpretation of the classical field is essential. 
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Several groups have recently observed the formation 
of a vortex lattice in weakly interacting Bose gases 
E3- El ■ Theoretically there have been attempts to 
understand the formation process 0,0,013 with simula- 
tions of the Gross-Pitaevskii equation for the condensate 
wavefunction. These papers either do not consider the 
effect of the noncondensed modes or they approxi- 
mate it with an added damping term. All of them stress 
the need for explicitly including this term since the vor- 
tices have to lose energy to be able to relax to a lattice 
configuration. 

Here we consider this problem within the framework 
of the classical theory of a complex field whose exact 
equation of motion is the standard nonlinear Schrodinger 
equation (NLSE). Within this framework we can give the 
true numerical answer which can be used as a testing 
ground for approximations such as the single mode one 
°f t° study the formation of the vortex lattice. If 
the dynamics of the field were ergodic we would expect 
that, in a fast enough rotating external potential, the 
vortex lattice would be formed without the need for any 
dissipative term. The issue of relaxation to thermal equi- 
librium for purely Hamiltonian dynamics is analogous to 
the Fermi-Ulam-Pasta problem |9( and has been studied 
also in the context of Bose gases 0, 0, 0, ^| . Here 
we study the formation of the lattice in three dimensions 
from an initially nonrotating Bose condensed gas both at 
T = and at finite temperatures using the NLSE. We 
find that, contrary to the prevalent view, the lattice is 
formed in both cases without any need for a damping 
term which suggests that thermalisation takes place in 
our system. 

We start our simulations with the nonrotating classi- 
cal field in thermal equilibrium. For T — the initial 
state is simply a pure condensate, that is, with a field 
proportional to the condensate wavefunction given by 
the Gross-Pitaevskii equation in the absence of rotation, 
ip = \/No(f> where iVo is the condensate atom number. 
For finite temperatures we sample the initial thermal 
equilibrium in the Bogoliubov approximation at a given 



temperature T for a fixed number iVo of condensate par- 
ticles. In this approximation the classical field is given by 
ip(r, 0) = \ZNq4>{t) + ip±(r). The random field tp±( T ) or- 
thogonal to <j) |l4j representing the thermal noise is given 
by 



(1) 



where the u n and v n are the Bogoliubov mode functions 
associated with <j> and the b n are independent random c- 
numbers taken from a Gaussian distribution that obeys 
the classical equipartition formula, {b^b n ) = kBT/e n , e n 
being the Bogoliubov energy of mode n. In practice, 
to sample this distribution we use the Brownian motion 
method described in |14| . In our work, the field tp is to be 
interpreted not as the condensate wavefunction but as the 
whole matter field (unlike 0000). We present here 
results from single realisations of the field ip which ex- 
perimentally correspond to single runs. We have checked 
that different realisations lead to similar results. 

In our simulations we consider a Bose condensed gas 
initially trapped in a cigar-shaped harmonic potential 
with oscillation frequencies whose ratio is 1:1:0.25, with 
10 5 atoms of mass m and a coupling constant g = 0.0343 
in units of huiaQ where u is the radial frequency and 
ao = y/ft/mu) is the oscillator length. The corresponding 
chemical potential is fi — 8huj. We start each simulation 
with the gas in thermal equilibrium. We abruptly turn 
on the trap anisotropy which leads to a change in the ra- 



dial frequencies: 



-x,y 



2 (l=Fe) where e = 0.025. Then 



the rotation frequency f2(f) of this anisotropy is slowly 
increased from zero to a final value £lf over 500W" 1 , to 
follow Procedure I in |l5f . After that we let the gas 
evolve in the presence of the rotating anisotropy until, 
for most of our simulations, the angular momentum of 
the gas reached a steady state. 

The calculation is performed in the rotating frame so 
that the NLSE takes the form 
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FIG. 1: Cut along the radial plane (z — 0) of the system 
spatial density at different times. Crosses (circles) indicate 
position of vortices of positive (negative) charge |17f . Left 
column: T = 0, fi/ = 0.8u. Top to bottom: initial state; near 
instability; turbulent behaviour; end of simulation. Right col- 
umn: ksT = 8huj, fi/ = 0.6lo. Top to bottom: initial state; 
entry of first vortex; entry of second vortex; end of simulation 
with a 3- vortex lattice. 
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FIG. 2: Angular momentum in units of fi per atom. The 
arrows marked E and C indicate the entry of the vortices 
into the condensate and the crystallisation of the lattice for 
fi/ = 0.8w. (a) T = 0, solid lines from bottom to top: fi//w = 
0.7(0), 0.75(7), 0.8(10); dashed line: fi//w = 0.8(10) with a 
gridsize of 64 x 64 x 256. All other curves were done on a 32 x 
32 x 128 grid. In parenthesis is the number of vortices in the 
lattice at the end of the simulation, (b) fesT = 4hui, flf/u) = 
0.4(0), 0.45(0), 0.5(0), 0.55(1), 0.6(1), 0.65(2), 0.7(6), 0.75(7), 
0.8(10). (c) k B T = 8hu>, Slf/u =0.4(0), 0.5(1), 0.55(1), 
0.6(3), 0.7(7), 0.75(7), 0.8(10). The arrows correspond to 
the approximate entry time of the vortices for fi/ = 0.6u> as 
shown in FigQ Note that the angular momentum shows no 
signature of the entries. 



where L z is the angular momentum operator along z and 
U is the anisotropic harmonic potential. The field ip is 
subject to periodic boundary conditions in the rotating 
frame. Our gridsize is 32 x 32 x 128 corresponding to an 
energy cutoff of 327iw per spatial direction, although we 
have also run simulations on a 64 x 64 x 256 grid (see be- 
low). We have checked that it is necessary to have time 



independent boundary conditions in the rotating frame: 
periodic boundary conditions in the lab frame arrest the 
rotation of the non-condensed gas. We have also checked 
that a pure condensate cannot be set into rotation by the 
effect of the boundary conditions only, since the conden- 
sate density is extremely weak at the borders of the grid 
for our choice of gridsize. The harmonic trap anisotropy 
is then a crucial element for the formation of the lattice 
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at T — by triggering the dynamic instability of . 

Zero initial temperature: This set of simulations can 
be divided into two groups: those for which the fi- 
nal rotation frequency is Vtf/uj < 0.7 and those with 
fi//w > 0.75. Between these two values lies the thresh- 
old for the dynamic instability of the condensate which 
changes the subsequent dynamics dramatically ^5|. In 
the first group, as the rotation frequency gradually in- 
creases with time, the condensate adiabatically follows a 
steady state, apart from excitations of the surface modes 
leading to a very small periodic variation of the angu- 
lar momentum (see curve for flf — 0.7a; in Fig|2Ji). 
With increasing ttf, the condensate's final state becomes 
more and more elliptically deformed, surrounded by a 
ring of vortices which however never enter it. The sec- 
ond group shows completely different behaviour when 
n(t)/uj ~ 0.75 (see left column of Fig2J: the instabil- 
ity sets in, the condensate becomes slightly S-shaped at 
t ~ 450a; -1 before being highly deformed and undergo- 
ing very turbulent motion [5j. This is accompanied by 
a large increase in angular momentum of the gas from 
almost zero when f2(t) < 0.75a; to between 5h-7ft, per 
particle (see FigEt) ■ At this point (t ~ 670a; -1 ) several 
vortices enter the high density region and, in less than 
200a; -1 , settle down to form a well-defined lattice. Af- 
ter this, a period of relaxation of around 800a; -1 begins 
with the initially rotating lattice finally stopping in the 
rotating frame. There remains a small random motion 
of the vortices around their equilibrium positions in the 
lattice together with density fluctuations in and around 
the condensate. 

At the end of the simulation, damping of the vortex 
motion has occurred and the initial energy of the vortex 
motion has been transferred in an effectively irreversible 
way to other degrees of freedom of the field. A similar 
phenomenon has been observed for the relative motion 
of two condensates [13j . If we assume that the field has 
reached a thermal distribution, we can calculate the tem- 
perature of the system by taking the final state of the 
simulation and evolving it with the conjugate gradient 
method in a trap rotating at f2/ . This reduces its energy 
and takes it to the local minimum associated with the 
vortex lattice. We then calculate the energy difference 
AE between the final state of the simulation and the one 
at the minimum. Assuming that Bogoliubov theory is 
valid, AE must correspond to the energy of a classical 
thermal distribution of weakly coupled harmonic oscilla- 
tors of amplitude b n which obeys the equipartition for- 
mula (6* b n )e n — k B T, with n being the Bogoliubov mode 
number. So, if Af is the number of modes in the system 
(and keeping in mind that we have to subtract the one 
corresponding to the condensate) then we have 

AE = 5>;& n )e„ = (AT - l)k B T. (3) 

n 

The final temperature is 0.6167iw for flf — 0.75a; and 



0.754?ia; for 17/ = 0.8a;, in other words it is extremely 
small, less than a tenth of the chemical potential. 

We have also carried out a simulation on a larger grid 
(64 x 64 x 256) to check the dependence on size. We 
chose ilf = 0.8a; and compared it with the one on the 
32 x 32 x 128 grid. The vortex nucleation and crystalli- 
sation phases are very similar and occur at roughly the 
same times. At longer times two differences arise: first, 
there are large underdamped oscillations of the angular 
momentum (see Fig(2K) • An analysis of the simulation 
suggests that these oscillations are those of the scissors 
mode. Second, the final temperature (0.094?ia;) differs by 
the ratio of the number of modes as expected: at time 
t = 500a; -1 when Vt(t) = Vtf, i/j had not yet reached the 
boundary in the smaller grid case and so the evolution 
of ip on both grids was identical up to this time with the 
same total energy which was conserved at later times re- 
sulting in the same value of AE. This exemplifies the fact 
that, in classical field theories, the relationship between 
energy and temperature depends on the energy cutoff. 

Since the thermal occupation of the modes is directly 
proportional to the temperature, we expect that all re- 
laxation processes which involve scattering from or into 
those modes (such as Landau-Beliaev damping) will be 
reduced. We are thus led to the conclusion that, for 
our simulations starting at T = 0, relaxation rates in 
the equilibration period after the formation of the lattice 
could depend on the size of the grid. However, with the 
present numerical results, we were not able to demon- 
strate this. 

Finite initial temperature: We performed simulations 
starting with fegT = Afiuj and fc^T = 8?kJ. Now not 
only the condensate but also other modes are occupied 
in the initial state, with a thermal distribution. For a 
final rotation frequency below that of the dynamic in- 
stability, the situation is quite different from that of the 
zero temperature case: the condensate is never deformed 
and the vortices do enter the condensate if flf > 0.55a; 
for k B T = Ahto and if % > 0.5a; for k B T = Shu. We 
have checked numerically that 0.51a; is the Landau crit- 
ical frequency above which the vortex-free condensate is 
no longer a local minimum of the energy. During the real 
time evolution of the right column of Figncorresponding 
to f2/ = 0.6a;, we find that the vortices enter only one at 
a time. That is, as the angular momentum of the cloud 
increases, one vortex, out of the group of vortices that 
surrounds the condensate will enter it and spiral slowly 
clockwise towards the center on a time scale of hundreds 
of a; -1 . After that vortex has reached the center, a sec- 
ond one enters slowly, repeating the trajectory of the first 
until it starts to interact with it and the two orbit around 
each other for a while after which a third will enter. At 
the end of the simulation, coinciding with the achieve- 
ment of the plateau in angular momentum, the lattice 
becomes stationary in the rotating frame and no further 
vortex enters the condensate (see right column of Fig^l. 
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For flf = 0.7uj we find that the condensate deforms itself 
elliptically after which three vortices enter at the same 
time and form a rotating lattice. After that, and spaced 
by several hundred w -1 , a fourth and then a fifth vortex 
enter. Finally, two further vortices enter simultaneously 
to form the final seven vortex lattice. At each interme- 
diate stage there is always a well defined lattice present 
although it is not stationary in the rotating frame. We 
should contrast this with the scenario of where a 
large number of vortices enter all at once into the conden- 
sate in a ring configuration and then some of them form 
a lattice while others are shed and leave the condensate. 

For flf above the dynamic instability frequency the 
situation is quite similar to the corresponding one at T = 
0. Once the instability has set in the lattice is formed for 
both temperatures in about 200lj _1 as in the T = case 
(see Fig|2]3,c-this weak temperature dependence was also 
found experimentally by |l8|). The main difference is 
that the relaxation time for the lattice to stop rotating is 
much shorter, on the order of a hundred not eight 
hundred. 

It is important to emphasize the multimode interpre- 
tation of the field. Transposing Penrose and Onsager's 
definition to the classical field theory, the condensate 
wavefunction is defined as the eigenvector corresponding 
to the largest eigenvalue of the one-body density matrix 
(^>*(r')^(r)) where the average is over an ensemble of 
initial states. If the system becomes turbulent because it 
encounters an instability, the trajectories of the neighbor- 
ing realisations will diverge exponentially. However, after 
averaging, we believe that the condensate wavefunction 
will not be a turbulent function. For T = there is only 
one initial state and so we replace ensemble averaging 
by one over time in the steady state regime |l9j . In our 
simulations with Clf = 0.8w the system must therefore 
be understood as becoming intrinsically multimode even 
though we started at T — with a pure condensate. 

Conclusions: We have identified two very different 
regimes for the crystallisation of the vortex lattice in the 
classical field theory. At T = when the dynamic insta- 
bility sets in (for high enough flf), the vortices enter the 
condensate abruptly and settle into a lattice even if there 
is no dissipative term in the equation of motion. This sce- 
nario is consistent with the one found experimentally at 
the ENS j2jj with comparable timescales. At finite tem- 
peratures, where 0,0 add an explicit damping term, we 
see that the thermal classical field can provide the damp- 
ing of the vortex motion on its own which leads to very 
different behaviour from that observed by those authors 
in that the vortices can enter one by one into the conden- 
sate and settle into a lattice before the entry of the follow- 
ing one. So far there has been no experimental check of 
this scenario. Finally, as the T — case shows, any theo- 
retical model which singles out the condensate mode for 



separate treatment with a Gross-Pitaevskii-type equation 
could run into trouble in turbulent situations since the 
separation between condensed and non-condensed modes 
would be hard to keep. 

We have been informed that crystallisation of the vor- 
tex lattice has also been observed in a simulation without 
a damping term by the group of N. Bigelow. We thank 
B. Durin, L. Carr, I. Carusotto, G. Shlyapnikov and J. 
Dalibard for useful contributions. C. L. acknowledges a 
fellowship from the Fundagao para a Ciencia e Tecnolo- 
gia of Portugal. LKB is a unit of ENS and of Universite 
Paris 6 associated to CNRS. 
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